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Abstract 

This paper introduces a novel algorithm to approximate the matrix with minimum nuclear 
norm among all matrices obeying a set of convex constraints. This problem may be understood as 
the convex relaxation of a rank minimization problem, and arises in many important applications 
as in the task of recovering a large matrix from a small subset of its entries (the famous Netflix 
problem). Off-the-shelf algorithms such as interior point methods are not directly amenable to 
large problems of this kind with over a million unknown entries. 

This paper develops a simple first-order and casy-to-implcmcnt algorithm that is extremely 
efficient at addressing problems in which the optimal solution has low rank. The algorithm is 
iterative and produces a sequence of matrices {X'^,Y''} and at each step, mainly performs a 
soft-thresholding operation on the singular values of the matrix Y'^. There are two remarkable 
features making this attractive for low-rank matrix completion problems. The first is that 
the soft-thresholding operation is applied to a sparse matrix; the second is that the rank of 
the iterates {X'^} is empirically nondecreasing. Both these facts allow the algorithm to make 
use of very minimal storage space and keep the computational cost of each iteration low. On 
the theoretical side, we provide a convergence analysis showing that the sequence of iterates 
converges. On the practical side, we provide numerical examples in which 1, 000 x 1, 000 matrices 
are recovered in less than a minute on a modest desktop computer. We also demonstrate that 
our approach is amenable to very large scale problems by recovering matrices of rank about 10 
with nearly a billion unknowns from just about 0.4% of their sampled entries. Our methods arc 
connected with the recent literature on linearized Bregman iterations for £i minimization, and 
we develop a framework in which one can understand these algorithms in terms of well-known 
Lagrange multiplier algorithms. 

Keywords. Nuclear norm minimization, matrix completion, singular value thresholding, La- 
grange dual function, Uzawa's algorithm. 



1 Introduction 

1.1 Motivation 

There is a rapidly growing interest in the recovery of an unknown low-rank or approximately low- 
rank matrix from very limited information. This problem occurs in many areas of engineering and 
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applied science such as machine learning [1,3,4], control [42] and computer vision, see [48]. As 
a motivating example, consider the problem of recovering a data matrix from a sampling of its 
entries. This routinely comes up whenever one collects partially filled out surveys, and one would 
like to infer the many missing entries. In the area of recommender systems, users submit ratings 
on a subset of entries in a database, and the vendor provides recommendations based on the user's 
preferences. Because users only rate a few items, one would like to infer their preference for unrated 
items; this is the famous Netflix problem [2]. Recovering a rectangular matrix from a sampling of 
its entries is known as the matrix completion problem. The issue is of course that this problem is 
extraordinarily ill posed since with fewer samples than entries, we have infinitely many completions. 
Therefore, it is apparently impossible to identify which of these candidate solutions is indeed the 
"correct" one without some additional information. 

In many instances, however, the matrix we wish to recover has low rank or approximately low 
rank. For instance, the Netflix data matrix of all user-ratings may be approximately low-rank 
because it is commonly believed that only a few factors contribute to anyone's taste or preference. 
In computer vision, inferring scene geometry and camera motion from a sequence of images is a well- 
studied problem known as the structure-from-motion problem. This is an ill-conditioned problem 
for objects may be distant with respect to their size, or especially for "missing data" which occur 
because of occlusion or tracking failures. However, when properly stacked and indexed, these images 
form a matrix which has very low rank (e.g. rank 3 under orthography) [21,48]. Other examples 
of low-rank matrix fitting abound; e.g. in control (system identification), machine learning (multi- 
class learning) and so on. Having said this, the premise that the unknown has (approximately) low 
rank radically changes the problem, making the search for solutions feasible since the lowest-rank 
solution now tends to be the right one. 

In a recent paper [13], Candes and Recht showed that matrix completion is not as ill-posed as 
people thought. Indeed, they proved that most low-rank matrices can be recovered exactly from 
most sets of sampled entries even though these sets have surprisingly small cardinality, and more 
importantly, they proved that this can be done by solving a simple convex optimization problem. 
To state their results, suppose to simplify that the unknown matrix M £ M"^" is square, and that 
one has available m sampled entries {Mij : (i,j) € il} where is a random subset of cardinality 
m. Then [13] proves that most matrices M of rank r can be perfectly recovered by solving the 
optimization problem 

minimize ll^ll* q -j^^, 

subject to Xij = Mij, {i,j) € 0., 

provided that the number of samples obeys 

m>Cn^/V log n (1.2) 

for some positive numerical constant In (1.1), the functional is the nuclear norm of the 

matrix M, which is the sum of its singular values. The optimization problem (1.1) is convex and 
can be recast as a semidefinite program [29,30]. In some sense, this is the tightest convex relaxation 
of the NP-hard rank minimization problem 

minimize rank(X) 
subject to Xij = Mij, {i,j) € ^, 

^Note that an n x n matrix of rank r depends upon r{2n — r) degrees of freedom. 
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since the nuclear ball {X : ||-X^||* < 1} is the convex hull of the set of rank-one matrices with 
spectral norm bounded by one. Another interpretation of Candes and Recht's result is that under 
suitable conditions, the rank minimization program (1.3) and the convex program (1.1) are formally 
equivalent in the sense that they have exactly the same unique solution. 

1.2 Algorithm outline 

Because minimizing the nuclear norm both provably recovers the lowest-rank matrix subject to 
constraints (see [45] for related results) and gives generally good empirical results in a variety of 
situations, it is understandably of great interest to develop numerical methods for solving (1.1). In 
[13], this optimization problem was solved using one of the most advanced semidefinite programming 
solvers, namely, SDPT3 [47]. This solver and others like SeDuMi are based on interior-point 
methods, and are problematic when the size of the matrix is large because they need to solve huge 
systems of linear equations to compute the Newton direction. In fact, SDPT3 can only handle 
n X n matrices with n < 100. Presumably, one could resort to iterative solvers such as the method 
of conjugate gradients to solve for the Newton step but this is problematic as well since it is well 
known that the condition number of the Newton system increases rapidly as one gets closer to the 
solution. In addition, none of these general purpose solvers use the fact that the solution may have 
low rank. We refer the reader to [40] for some recent progress on interior-point methods concerning 
some special nuclear norm-minimization problems. 

This paper develops the singular value thresholding algorithm for approximately solving the 
nuclear norm minimization problem (1.1) and by extension, problems of the form 

minimize ll^ll* 
subject to -4(^) = 

where ^ is a linear operator acting on the space of ni x 712 matrices and b £ M™. This algorithm is 
a simple first-order method, and is especially well suited for problems of very large sizes in which 
the solution has low rank. We sketch this algorithm in the special matrix completion setting and 
let Vn be the orthogonal projector onto the span of matrices vanishing outside of Q so that the 
(i,j)th component of 'Pn(X) is equal to Xij if (i,j) G and zero otherwise. Our problem may be 
expressed as 

minimize ||-'^||* /-, t;^ 

subject to Vn{X)=Vn{M), ^ 

with optimization variable X G ]K"iX"2_ pj^ 7- > g and a sequence {5k}k>i of scalar step sizes. 
Then starting with = G R'"!^"^^ the algorithm inductively defines 

X'' = shrink(l"''~\r), 

= Y''-^ + 6kVniM - X'') 

until a stopping criterion is reached. In (1.6), shrink(l^,r) is a nonlinear function which applies a 
soft-thresholding rule at level r to the singular values of the input matrix, see Section 2 for details. 
The key property here is that for large values of r, the sequence {X^} converges to a solution which 
very nearly minimizes (1.5). Hence, at each step, one only needs to compute at most one singular 
value decomposition and perform a few elementary matrix additions. Two important remarks are 
in order: 
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1. Sparsity. For each k > 0, vanishes outside of Q. and is, therefore, sparse, a fact which can 
be used to evaluate the shrink function rapidly. 

2. Low-rank property. The matrices turn out to have low rank, and hence the algorithm has 
minimum storage requirement since we only need to keep principal factors in memory. 

Our numerical experiments demonstrate that the proposed algorithm can solve problems, in 
Matlab, involving matrices of size 30, 000 x 30, 000 having close to a billion unknowns in 17 minutes 
on a standard desktop computer with a 1.86 GHz CPU (dual core with Matlab's multithreading 
option enabled) and 3 GB of memory. As a consequence, the singular value thresholding algorithm 
may become a rather powerful computational tool for large scale matrix completion. 

1.3 General formulation 

The singular value thresholding algorithm can be adapted to deal with other types of convex 
constraints. For instance, it may address problems of the form 

minimize ||^||* q 
subject to fi{^) < 0, i = 1, . . . , m, 

where each fi is a Lipschitz convex function (note that one can handle linear equality constraints by 
considering pairs of affine functionals). In the simpler case where the /j's are affine functionals, the 
general algorithm goes through a sequence of iterations which greatly resemble (1.6). This is useful 
because this enables the development of numerical algorithms which are effective for recovering 
matrices from a small subset of sampled entries possibly contaminated with noise. 

1.4 Contents and notations 

The rest of the paper is organized as follows. In Section 2, we derive the singular value threshold- 
ing (SVT) algorithm for the matrix completion problem, and recasts it in terms of a well-known 
Lagrange multiplier algorithm. In Section 3, we extend the SVT algorithm and formulate a gen- 
eral iteration which is applicable to general convex constraints. In Section 4, we establish the 
convergence results for the iterations given in Sections 2 and 3. We demonstrate the performance 
and effectiveness of the algorithm through numerical examples in Section 5, and review additional 
implementation details. Finally, we conclude the paper with a short discussion in Section 6. 

Before continuing, we provide here a brief summary of the notations used throughout the 
paper. Matrices are bold capital, vectors are bold lowercase and scalars or entries are not bold. 
For instance, X is a matrix and Xij its (i,j)th entry. Likewise, x is a vector and Xi its ith 
component. The nuclear norm of a matrix is denoted by ||X||*, the Frobenius norm by 
and the spectral norm by ||-X^||2; note that these are respectively the 1-norm, the 2-norm and the 
sup-norm of the vector of singular values. The adjoint of a matrix X is X* and similarly for 
vectors. The notation diag(a;), where a; is a vector, stands for the diagonal matrix with {xi} as 
diagonal elements. We denote by {X,Y) = trace(X*l^) the standard inner product between two 
matrices (||X|||^ = {X,X)). The Cauchy-Schwarz inequality gives {X,Y) < \\X\\f\\Y\\p and it is 
well known that we also have {X,Y) < || X ||*||1^||2 (the spectral and nuclear norms are dual from 
one another), see e.g. [13,45]. 
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2 The Singular Value Thresholding Algorithm 



This section introduces the singular value thresholding algorithm and discusses some of its ba- 
sic properties. We begin with the definition of a key building block, namely, the singular value 
thresholding operator. 

2.1 The singular value shrinkage operator 

Consider the singular value decomposition (SVD) of a matrix X G ]R"i^"2 gf rank r 

X = UT.V\ S = diag({(7,}i<,<,), (2.1) 

where U and V are respectively ni x r and n2 x r matrices with orthonormal columns, and the 
singular values dj are positive (unless specified otherwise, we will always assume that the SVD of 
a matrix is given in the reduced form above). For each r > 0, we introduce the soft-thresholding 
operator defined as follows: 

Vr{X) := UVr{^)V*, Vr{T.) = diag({c7i - r)+}), (2.2) 

where t+ is the positive part of t, namely, t+ = max(0, t). In words, this operator simply applies a 
soft-thresholding rule to the singular values of X, effectively shrinking these towards zero. This is 
the reason why we will also refer to this transformation as the singular value shrinkage operator. 
Even though the SVD may not be unique, it is easy to see that the singular value shrinkage operator 
is well defined and we do not elaborate further on this issue. In some sense, this shrinkage operator 
is a straightforward extension of the soft-thresholding rule for scalars and vectors. In particular, 
note that if many of the singular values of X are below the threshold r, the rank of T>-r{X) may 
be considerably lower than that of X, just like the soft-thresholding rule applied to vectors leads 
to sparser outputs whenever some entries of the input are below threshold. 

The singular value thresholding operator is the proximity operator associated with the nuclear 
norm. Details about the proximity operator can be found in e.g. [35]. 

Theorem 2.1 For each r > and Y £ ]K"iX"2^ ^/jg singular value shrinkage operator (2.2) obeys 

Vr{Y) = a.Tgmm !^^\\X -Y\\l + T\\X\Uy (2.3) 

Proof. Since the function hQ{X) := t||X||^, -|- ^\\X — YWj^ is strictly convex, it is easy to see that 
there exists a unique minimizer, and we thus need to prove that it is equal to T>t-{Y). To do this, 
recall the definition of a subgradient of a convex function / : M"i^"2 ]R. We say that Z is a 
subgradient of / at Xq, denoted Z £ 5/(Xo), if 

f{X)>f{Xo) + {Z,X-Xo) (2.4) 

for all X. Now X minimizes Hq if and only if is a subgradient of the functional Hq at the point 
X, i.e. 

0€ X -Y + Td\\X\U, (2.5) 

where is the set of subgradients of the nuclear norm. Let X G ]R'^iX"2 |-,g arbitrary 

matrix and VSV* be its SVD. It is known [13,37,49] that 

d\\X\\, = {UV* + W : W U*W = 0, WV = 0, ||W||2 < l} • (2.6) 
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Set X := Dt-IY) for short. In order to show that X obeys (2.5), decompose the SVD of Y as 



Y = C/oSo^o* + Ui^iV{, 

where Uq, Vq (resp. Ui, Vi) are the singular vectors associated with singular values greater than r 
(resp. smaller than or equal to r). With these notations, we have 

X = Uoi^o - rI)Vo* 

and, therefore, 

Y - X = t{UoVq* + W), W = T-^Uii:iV{ . 

By definition, UqW = 0, WVq = and since the diagonal elements of Si have magnitudes 
bounded by r, we also have ||VF||2 < 1. Hence Y — X € Td\\X\\^:, which concludes the proof. q 

2.2 Shrinkage iterations 

We are now in the position to introduce the singular value thresholding algorithm. Fix r > and 
a sequence {6^} of positive step sizes. Starting with Yq, inductively define for /c = 1, 2, . . ., 

ix^ = vAY''-^), 

l^fe = r'^-i + 6kVn{M -X^) 

until a stopping criterion is reached (we postpone the discussion this stopping criterion and of the 
choice of step sizes). This shrinkage iteration is very simple to implement. At each step, we only 
need to compute an SVD and perform elementary matrix operations. With the help of a standard 
numerical linear algebra package, the whole algorithm can be coded in just a few lines. 

Before addressing further computational issues, we would like to make explicit the relationship 
between this iteration and the original problem (1.1). In Section 4, we will show that the sequence 
{X^} converges to the unique solution of an optimization problem closely related to (1.1), namely. 



minimize t||X||* + , 
subject to Vn{X)=Vn{M). ^ ' 



Furthermore, it is intuitive that the solution to this modified problem converges to that of (1.5) as 
r — > oo as shown in Section 3. Thus by selecting a large value of the parameter r, the sequence of 
iterates converges to a matrix which nearly minimizes (1.1). 

As mentioned earlier, there are two crucial properties which make this algorithm ideally suited 
for matrix completion. 

• Low-rank property. A remarkable empirical fact is that the matrices in the sequence {X^} 
have low rank (provided, of course, that the solution to (2.8) has low rank). We use the word 
"empirical" because all of our numerical experiments have produced low-rank sequences but 
we cannot rigorously prove that this is true in general. The reason for this phenomenon is, 
however, simple: because we are interested in large values of r (as to better approximate the 
solution to (1.1)), the thresholding step happens to 'kill' most of the small singular values 
and produces a low-rank output. In fact, our numerical results show that the rank of X^ is 
nondecreasing with k, and the maximum rank is reached in the last steps of the algorithm, 
see Section 5. 



6 



Thus, when the rank of the solution is substantiahy smaher than either dimension of the 
matrix, the storage requirement is low since we could store each in its SVD form (note 
that we only need to keep the current iterate and may discard earlier values). 

• Sparsity. Another important property of the SVT algorithm is that the iteration matrix 
is sparse. Since 1^*^ = 0, we have by induction that vanishes outside of Q. The fewer 
entries available, the sparser Y''. Because the sparsity pattern Q is fixed throughout, one can 
then apply sparse matrix techniques to save storage. Also, if \^}\ = m, the computational cost 
of updating Y^ is of order m. Moreover, we can call subroutines supporting sparse matrix 
computations, which can further reduce computational costs. 

One such subroutine is the SVD. However, note that we do not need to compute the entire 
SVD of Y^ to apply the singular value thresholding operator. Only the part corresponding 
to singular values greater than r is needed. Hence, a good strategy is to apply the iterative 
Lanczos algorithm to compute the first few singular values and singular vectors. Because 
Y^ is sparse, Y'' can be applied to arbitrary vectors rapidly, and this procedure offers a 
considerable speedup over naive methods. 

2.3 Relation with other works 

Our algorithm is inspired by recent work in the area of ii minimization, and especially by the work 
on linearized Bregman iterations for compressed sensing, see [9-11,23,44,51] for linearized Bregman 
iterations and [14-17, 26] for some information about the field of compressed sensing. In this line 
of work, linearized Bregman iterations are used to find the solution to an underdetermined system 
of linear equations with minimum ii norm. In fact, Theorem 2.1 asserts that the singular value 
thresholding algorithm can be formulated as a linearized Bregman iteration. Bregman iterations 
were first introduced in [43] as a convenient tool for solving computational problems in the imaging 
sciences, and a later paper [51] showed that they were useful for solving ^i-norm minimization 
problems in the area of compressed sensing. Linearized Bregman iterations were proposed in [23] 
to improve performance of plain Bregman iterations, see also [51]. Additional details together 
with a technique for improving the speed of convergence called kicking are described in [44]. On 
the practical side, the paper [11] applied Bregman iterations to solve a deblurring problem while 
on the theoretical side, the references [9, 10] gave a rigorous analysis of the convergence of such 
iterations. New developments keep on coming out at a rapid pace and recently, [32] introduced a 
new iteration, the split Bregman iteration, to extend Bregman-type iterations (such as linearized 
Bregman iterations) to problems involving the minimization of -^i-like functionals such as total- 
variation norms, Besov norms, and so forth. 

When applied to ^i-minimization problems, linearized Bregman iterations are sequences of 
soft-thresholding rules operating on vectors. Iterative soft-thresholding algorithms in connection 
with ii or total-variation minimization have quite a bit of history in signal and image processing 
and we would like to mention the works [12, 39] for total-variation minimization, [24, 25, 31] for 
li minimization, and [5, 7, 8, 19, 20, 27, 28, 46] for some recent applications in the area of image 
inpainting and image restoration. Just as iterative soft-thresholding methods are designed to find 
sparse solutions, our iterative singular value thresholding scheme is designed to find a sparse vector 
of singular values. In classical problems arising in the areas of compressed sensing, and signal or 
image processing, the sparsity is expressed in a known transformed domain and soft-thresholding is 
applied to transformed coefficients. In contrast, the shrinkage operator Vr is adaptive. The SVT not 
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only discovers a sparse singular vector but also the bases in which we have a sparse representation. 
In this sense, the SVT algorithm is an extension of earlier iterative soft-thresholding schemes. 

Finally, we would like to contrast the SVT iteration (2.7) with the popular iterative soft- 
thresholding algorithm used in many papers in imaging processing and perhaps best known under 
the name of Proximal Forward-Backward Splitting method (PFBS), see [8,22,24,31,33] for example. 
The constrained minimization problem (1.5) may be relaxed into 



for some A > 0. Theorem 2.1 asserts that Dx is the proximity operator of A||X||* and Proposition 
3.1(iii) in [22] gives that the solution to this unconstrained problem is characterized by the fixed 
point equation X = Dxsi^ + SPq{M — X)) for each 5 > 0. One can then apply a simplified 
version of the PFBS method (see (3.6) in [22]) to obtain iterations of the form 



The difference with (2.7) may seem subtle at first — replacing X'' in (2.10) with "J^^-i and setting 
= ^ gives (2.7) with r = \5 — but has enormous consequences as this gives entirely different 
algorithms. First, they have different limits: while (2.7) converges to the solution of the constrained 
minimization (2.8), (2.10) converges to the solution of (2.9) provided that the sequence of step sizes 
is appropriately selected. Second, selecting a large A (or a large value of r = X6) in (2.10) gives 
a low-rank sequence of iterates and a limit with small nuclear norm. The limit, however, does 
not fit the data and this is why one has to choose a small or moderate value of A (or of r = X6). 
However, when A is not sufficiently large, the X^ may not have low rank even though the solution 
has low rank (and one may need to compute many singular vectors), and is not sufficiently 
sparse to make the algorithm computationally attractive. Moreover, the limit does not necessary 
have a small nuclear norm. These are reasons why (2.10) is not suitable for matrix completion. 

2.4 Interpretation as a Lagrange multiplier method 

In this section, we recast the SVT algorithm as a type of Lagrange multiplier algorithm known as 
Uzawa's algorithm. An important consequence is that this will allow us to extend the SVT algorithm 
to other problems involving the minimization of the nuclear norm under convex constraints, see 
Section 3. Further, another contribution of this paper is that this framework actually recasts linear 
Bregman iterations as a very special form of Uzawa's algorithm, hence providing fresh and clear 
insights about these iterations. 

In what follows, we set fr{X) = t\\X\U + ^\\X\\l for some fixed r > and recall that we wish 



minimize A||X||, + -\\Vn{X) - Vn{M)fp 



(2.9) 



Introducing an intermediate matrix Y , this algorithm may be expressed as 




(2.10) 



to solve (2.8) 



minimize /r(-X^) 

subject to Vn{X) = Vn{M). 



The Lagrangian for this problem is given by 



CiX,Y 



) = U{X) + {Y,Vn{M-X)) 
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where Y E ]R"i^"2_ Strong duality holds and X* and Y* are primal-dual optimal if {X* , Y*) is a 
saddlepoint of the Lagrangian C{X, Y), i.e. a pair obeying 

supinf£(X,r) =C{X*,Y'') = inf sup£(X, r). (2.11) 

(The function go{Y) = mix ,Y) is called the dual function.) Uzawa's algorithm approaches 
the problem of finding a saddlepoint with an iterative procedure. From Yo = 0, say, inductively 
define 

f £(X^ l^'^-i) = minx JC{X, Y^-^) 
\y>' = Y>'-^ + 6krn{M -Xf"), 

where {Sk}k>i is a sequence of positive step sizes. Uzawa's algorithm is, in fact, a subgradient 
method applied to the dual problem, where each step moves the current iterate in the direction of 
the gradient or of a subgradient. Indeed, observe that 

dY9o{Y) = dYC{X,Y)=Vn{M-X), (2.13) 

where X is the minimizer of the Lagrangian for that value of Y so that a gradient descent update 
for Y is of the form 

yk = + 5kdY9o{Y''~^) = r'^-i + 6kVn{M - X''). 
It remains to compute the minimizer of the Lagrangian (2.12), and note that 

argmin /,(X) + {Y,VniM - X)) = argmin r||X||, + ^\\X - VnYfp. (2.14) 

However, we know that the minimizer is given by 'DT-{T'fi(Y)) and since Y'' = Vn{Y'') iov all k > 0, 
Uzawa's algorithm takes the form 

jx'' =VriY''-'^) 

which is exactly the update (2.7). This point of view brings to bear many different mathemat- 
ical tools for proving the convergence of the singular value thresholding iterations. For an early 
use of Uzawa's algorithm minimizing an ^i-like functional, the total-variation norm, under linear 
inequality constraints, see [12]. 

3 General Formulation 

This section presents a general formulation of the SVT algorithm for approximately minimizing the 
nuclear norm of a matrix under convex constraints. 

3.1 Linear equality constraints 

Set the objective functional /r(X) = ''"ll-X^II* + ^||^|||' for some fixed r > 0, and consider the 
following optimization problem: 

minimize friX) 

subject to A{X) = b, ^ ' ^ 



9 



where ^ is a linear transformation mapping rii x n2 matrices into ffi"* („4* is the adjoint of A). This 
more general formulation is considered in [13] and [45] as an extension of the matrix completion 
problem. Then the Lagrangian for this problem is of the form 

C{X,y) = U{X) + {y,b-A{X)), (3.2) 

where X £ J^'^ix'^a g^j^^^ y £ j^m^ g^j^^^^ starting with = 0, Uzawa's iteration is given by 

yk-^+6k{b-A{X'^)). 

The iteration (3.3) is of course the same as (2.7) in the case where ^ is a sampling operator 
extracting m entries with indices in Q out of an rii x n2 matrix. To verify this claim, observe 
that in this situation, A* A = Vn, and let M be any matrix obeying A{M) = b. Then defining 
= A*{y'') and substituting this expression in (3.3) gives (2.7). 

3.2 General convex constraints 

One can also adapt the algorithm to handle general convex constraints. Suppose we wish to 
minimize /r(^) defined as before over a convex set X £ C. To simplify, we will assume that this 
convex set is given by 

C = {X : n{X) < 0, Vi = l,...,m}, 

where the /j's are convex functionals (note that one can handle linear equality constraints by 
considering pairs of affine functionals). The problem of interest is then of the form 

minimize fr{X) ,^ 
subject to fi{^) < 0, i = 1, . . . , m. 

Just as before, it is intuitive that as r — > oo, the solution to this problem converges to a minimizer 
of the nuclear norm under the same constraints (1.7) as shown in Theorem 3.1 at the end of this 
section. 

Put J-{X) := (/i(X), . . . , fm{X)) for short. Then the Lagrangian for (3.4) is equal to 

£iX,y)=fr{X) + {y,J^{X)), 

where X £ ]^"i><"2 gj^^j y ^ j^m -g ^ vector with nonnegative components denoted, as usual, 
by y > 0. One can apply Uzawa's method just as before with the only modification that we will 
use a subgradient method with projection to maximize the dual function since we need to make 
sure that the successive updates y'^ belong to the nonnegative orthant. This gives 

(X'^ = argmin{^(X) + {y''~\ J^{X))}, 
\y>^ = [y>'~^+6kHX'')] + - 

Above, a;_|_ is of course the vector with entries equal to max(xi,0). When is an affine mapping 
of the form b — A{X) so that one solves 

minimize fri^) 
subject to A{X) > b, 
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this simplifies to 

[y^~^+5k{b-A{X^))]+, 
and thus the extension to hnear inequahty constraints is straightforward. 

3.3 Example 

An interesting example concerns the extension of the Dantzig selector [18] to matrix problems. 
Suppose we have available linear measurements about a matrix M of interest 

h = A{M)+z, (3.7) 

where z G R'" is a noise vector. Then under these circumstances, one might want to find the 
matrix which minimizes the nuclear norm among all matrices which are consistent with the data b. 
Inspired by the work on the Dantzig selector which was originally developed for estimating sparse 
parameter vectors from noisy data, one could approach this problem by solving 

minimize ||^||* /o q\ 

subject to |vec(^*(r))| < wec{E), r := b - AiX), ^ ' 

where E is an array of tolerances, which is adjusted to fit the noise statistics [18]. Above, vec(A) < 
vec(S), for any two matrices A and B, means componentwise inequalities; that is, Aij < Bij for all 
indices i, j. We use this notation as not to confuse the reader with the positive semidefinite ordering. 
In the case of the matrix completion problem where A extracts sampled entries indexed by 0, one 
can always see the data vector as the sampled entries of some matrix B obeying V^{B) = A*{b); 
the constraint is then natural for it may be expressed as 

\Bij - Xij\ < Eij, G n, 

If z is white noise with standard deviation a, one may want to use a multiple of a for Eij. In 
words, we are looking for a matrix with minimum nuclear norm under the constraint that all of its 
sampled entries do not deviate too much from what has been observed. 

Let y+ G i^nixn2 (resp. Y- £ i^"-iX"-2^ ^j^g Lagrange multiplier associated with the compo- 
nentwise linear inequality constraints \'ec{A*{r)) < vec{E) (resp. —vec{A*{r)) < vec{E)). Then 
starting with = 0, the SVT iteration for this problem is of the form 



fx^' = Vr{A*AiY^~' - Yi-i)), 

[Y^ = [Yi~^ + 6k{±A*ir'') - E)]+, r*^ = 6'^ - AiX''), ^ ' ' 

where again [•]+ is applied componentwise. 

We conclude by noting that in the matrix completion problem where A* A = Vfi and one 
observes V^{B), one can check that this iteration simplifies to 



Vr{Y^'^ -Y'I-^), 

[Yt^ + 6kVni±iB - X'^) - E)]+. 



Again, this is easy to implement and whenever the solution has low rank, the iterates X^ have low 
rank as well. 
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3.4 When the proximal problem gets close 

We now show that minimizing the proximal objective /r(^) = tH-X^H* + is the same as 

minimizing the nuclear norm in the limit of large r's. The theorem below is general and covers the 
special case of linear equality constraints as in (2.8). 

Theorem 3.1 Let X* be the solution to (3.4) and X^q be the minimum Frobenius-norm solution 
to (1.7) defined as 

Xq^ : = arg min is a solution of (1.7)}. (3-11) 

Assume that the fi{X) 's, 1 < i < m, are convex and lower semi- continuous. Then 

lim ||X* - XooIIf = 0. (3.12) 

T— >00 

Proof. It follows from the definition of X* and X^o that 

II-^tII* ~^ 7;~II-^tIIf — l|-^oo||* + 7;~ ll-^oo IIf; &T^d \\Xoo\\* ^ ll-^rll*- (3.13) 
zr 2t 

Summing these two inequalities gives 

< ii^ooIIf, (3.14) 

which implies that ||J5C*|||, is bounded uniformly in r. Thus, we would prove the theorem if we 
could establish that any convergent subsequence {X*^}k>i must converge to X^q. 

Consider an arbitrary converging subsequence {X*^} and set Xc ■= limi^^^ X*^. Since for 
each 1 < i < m, fi{X*^) < and fi is lower semi-continuous, Xc obeys 

fiiXc)<0, i = l,...,m. (3.15) 

Furthermore, since \\X*fp is bounded, (3.13) yields 

limsup < ||Xoo||*, ll^ooll* < liminf 

T^oo r^oo 

An immediate consequence is limT-_»oo ||^*||* = ||-X^oo||* and, therefore, ||^c||* = ||^oo||*- This 
shows that Xc is a solution to (1.1). Now it follows from the definition of X^o that ||JCc||f > 
||-X^oo||f, while we also have ||Xc||_f < ||Xoo||f because of (3.14). We conclude that ||Xc||f = 
||Xoo||f and thus Xc = Xqo since Xoo is unique. □ 



4 Convergence Analysis 

This section establishes the convergence of the SVT iterations. We begin with the simpler proof 
of the convergence of (2.7) in the special case of the matrix completion problem, and then present 
the argument for the more general constraints (3.5). We hope that this progression will make the 
second and more general proof more transparent. 
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4.1 Convergence for matrix completion 

We begin by recording a lemma which estabhshes the strong convexity of the objective fr- 

Lemma 4.1 Let Z £ dfr{X) and Z' £ dfr{X'). Then 

{Z - Z\X - X') >\\X - X'fp. (4.1) 

Proof. An element Z of dfriX) is of the form Z = tZq + X, where Zq G (9|| and similarly 
for Z' . This gives 

{Z -Z' X- X') = t{Zo- Z'q, X-X') + \\X - X'fp 

and it thus suffices to show that the first term of the right-hand side is nonnegative. From (2.6), 
we have that any subgradient of the nuclear norm at X obeys ||-^o||2 ^ 1 and {Zq, X) = ||^||*. In 
particular, this gives 

\{Zq,X')\ < ||Zo||2||^'||* < ll-X^'ll^,, \{Zq,X)\ < \\Zq\\2\\X\\^ < \\x\\^. 

Whence, 

{Zo -Z'o,X- X') = (Zo, X) + {Zi X') - (Zo, X') - {Z'o, X) 
= \\x\u + \\x'\u - {Zo,X') - {Z^X) > 0, 

which proves the lemma. □ 
This lemma is key in showing that the SVT algorithm (2.7) converges. 

Theorem 4.2 Suppose that the sequence of step sizes obeys < inf (5fc < sup^fc < 2. Then the 
sequence {X^} obtained via (2.7) converges to the unique solution of (2.8). 

Proof. Let {X*,Y*) be primal-dual optimal for the problem (2.8). The optimality conditions give 

= Z'' - Vn{Y''~'^) 
= Z" -VniY*), 

for some G dfr{X^) and some Z* £ dfr{X*). We then deduce that 

[Z^ - Z*) - Vn{Y^^^ -Y*) = Q 

and, therefore, it follows from Lemma 4.1 that 

{X^ - X*, Vn{Y^^^ - Y^)) = (Z^ - Z*, X^ - X^) > \\X^ - X*fp. (4.2) 

We continue and observe that because VnX* = VnM, 

WVniY'^ - Y^)\\f = WVniY'^-' - Y*) + 6kVn{X* - X^)\\f. 

Therefore, setting ru = \\Vn{Y'' - 1"*)||f, 

rl = rl_, - 25k{Vn{Y^-^ - Y*),X^ - X*) + 5l\\Vn{X* - X^)\\l 

< rl_^ - 25k\\X^ - X*fp + 5l\\X^ - X*fp (4.3) 

since for any matrix X, \\Vn{X)\\F < Under our assumptions about the size of 5^, we have 

26k — 5\> [3 for all A; > 1 and some /5 > and thus 

rl<rl^,-P\\X'-X^fp. (4.4) 

Two properties follow from this: 
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1. The sequence 

2. As a consequence, WX'^ 
The theorem is estabhshed. 

4.2 General convergence theorem 

Our second result is more general and establishes the convergence of the SVT iterations to the 
solution of (3.4) under general convex constraints. From now now, we will only assume that the 
function J-{X) is Lipschitz in the sense that 

<L(^)||X-l^||i., (4.5) 

for some nonnegative constant L{T). Note that if T is afhne, J'{X) = b — A{X), we have 
L(JF) = 1 1 2 where ||.4||2 is the spectrum norm of the linear transformation A defined as ||.4||2 := 
sup{||^(X)||£2 : = !}• We also recall that J^iX) = {fi{X), . . . , fm{X)) where each fi is 

convex, and that the Lagrangian for the problem (3.4) is given by 

C{X,y) = fr{X) + {y,J^{X)), y > 0. 

We will assume to simplify that strong duality holds which is automatically true if the constraints 
obey constraint qualifications such as Slater's condition [6]. 
We first establish the following preparatory lemma. 

Lemma 4.3 Let {X*,y*) be a primal-dual optimal pair for (3.4). Then for each 5 > 0, y* obeys 

y* = [y* + 6nx*)U. (4.6) 
Proof. Recall that the projection xq of a point x onto a convex set C is characterized by 

(xo G C, 

\{y - xo,x - xq) <0, yy € C. 
In the case where C = = {tc G : x > 0}, this condition becomes Xq > and 

{y - xo,x- Xq) < 0, Vy > 0. 
Now because y* is dual optimal we have 

C{X\y*)>C{X\y), Vy > 0. 
Substituting the expression for the Lagrangian, this is equivalent to 

{y-y\T{X*)) <0, Vy>0, 

which is the same as 

(y - y^ y* + pT{X^) - y*) < 0, Vy > 0, Vp > 0. 

Hence it follows that y* must be the projection of y* + pJ-{X*) onto the nonnegative orthant W^. 
Since the projection of an arbitrary vector x onto is given by 0^+, our claim follows. q 

We are now in the position to state our general convergence result. 



l^*)||i?} is nonincreasing and, therefore, converges to a limit. 
^ as A; ^ oo. 



□ 
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Theorem 4.4 Suppose that the sequence of step sizes obeys < inf 5^ < supj^ < 2/||L(J^)p, 
where L{T) is the Lipschitz constant in (4.5). Then assuming strong duality, the sequence {-X^*^} 
obtained via (3.5) converges to the unique solution of (3.4). 

Proof. Let {X*,y*) be primal-dual optimal for the problem (3.4). We claim that the optimality 
conditions give that for all X 

{Z'' ,X - X^) + {y''-\j^{X) - F{X^)) >0, 

{Z\ X-X*) + {y\r{X) - J^iX")) > 0, (4.7) 

for some Z'' G dfr{X^) and some Z* G dfr{X*). We justify this assertion by proving one of the 
two inequalities since the other is exactly similar. For the first, X'' minimizes vC(JC, j/'^"^) over all 
X and, therefore, there exist Z'' G dfr{X^) and Z^ G dfi{X^), 1 < i < m, such that 



1=1 

Now because each fi is convex, 

/,(X)-/,(X'=)>(Zf,X-X'=) 

and, therefore. 



{z\x -x'') + J2 yt\Mx) - f.ix'^)) > + E yt'zt X - x') = o. 

i=l i=l 

This is (4.7). 

Now write the first inequality in (4.7) for X*, the second for X'^ and sum the two inequalities. 
This gives 

(Z^ - Z\ X^ - X*) + (/"I - y*, T{X^) - T{X*)) < 0. 

The rest of the proof is essentially the same as that of Theorem 4.5. It follows from Lemma 4.1 
that 

{y'"^^ - y\Hx'') - Hx*)) < -{z^ - z\x'' -x*) < -\\x^ -x*\\l. (4.8) 

We continue and observe that because y* = [y* + 6k^{X)]^ by Lemma 4.3, we have 

Wy'^ -y*\\ = \\[y^-' + dfc.^(X^)]+ - [y* + <5fc.F(X^)]+|| 
< \\y'~'-y* + Sk{HX'')-HXm 

since the projection onto the convex set is a contraction. Therefore, 

Wy" - y*f = Wy"-' - y*f + 2Sk {y""' - y\T{x^) - Hx")) + 5lll^(x^') - Hx'')f 

< 11/"^ - y*f - 25fc||^'' - X*fp + 5lL^ \\x^ - 

where we have put L instead of L{T) for short. Under our assumptions about the size of 5^, we 
have 25k - ^II"^ > /3 for ah fc > 1 and some /3 > 0. Then 

11/ - 2/11' < Wy"-' - y*f - Plix"^ - x^Wl, (4.9) 
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and the conclusion is as before. 

The problem (3.1) with linear constraints can be reduced to (3.4) by choosing 



" b 




' A ' 


-b 




-A 



and we have the following corollary: 

Corollary 4.5 Suppose that the sequence of step sizes obeys < ini6k < supJ^ < 2/||^||2. Then 
the sequence {X''} obtained via (3.3) converges to the unique solution of (3.1). 

Let PII2 := sup{P(X)||f : \\X\\f = 1}. With J^(X) given as above, we have \L{T)\^ = 2\\A\\l 
and thus, Theorem 4.4 guarantees convergence as long as < iniSk < supiJ^ < 1/||^||2- However, 
an argument identical to the proof of Theorem 4.2 would remove the extra factor of two. We omit 
the details. 



5 Implementation and Numerical Results 

This section provides implementation details of the SVT algorithm — as to make it practically 
effective for matrix completion — such as the numerical evaluation of the singular value thresholding 
operator, the selection of the step size 5k, the selection of a stopping criterion, and so on. This 
section also introduces several numerical simulation results which demonstrate the performance 
and effectiveness of the SVT algorithm. We show that 30, 000 x 30, 000 matrices of rank 10 are 
recovered from just about 0.4% of their sampled entries in a matter of a few minutes on a modest 
desktop computer with a 1.86 GHz CPU (dual core with Matlab's multithreading option enabled) 
and 3 GB of memory. 

5.1 Implementation details 

5.1.1 Evaluation of the singular value thresholding operator 

To apply the singular value tresholding operator at level r to an input matrix, it suffices to know 
those singular values and corresponding singular vectors above the threshold r. In the matrix 
completion problem, the singular value thresholding operator is applied to sparse matrices {^^'^j 
since the number of sampled entries is typically much lower than the number of entries in the 
unknown matrix M, and we are hence interested in numerical methods for computing the dominant 
singular values and singular vectors of large sparse matrices. The development of such methods is 
a relatively mature area in scientific computing and numerical linear algebra in particular. In fact, 
many high-quality packages are readily available. Our implementation uses PROPACK, see [36] 
for documentation and availability. One reason for this choice is convenience: PROPACK comes 
in a Matlab and a Fortran version, and we find it convenient to use the well-documented Matlab 
version. More importantly, PROPACK uses the iterative Lanczos algorithm to compute the singular 
values and singular vectors directly, by using the Lanczos bidiagonalization algorithm with partial 
reorthogonalization. In particular, PROPACK does not compute the eigenvalues and eigenvectors 
of {y'')*Y'' and Y^{Y'')*, or of an augmented matrix as in the Matlab built-in function 'svds' for 
example. Consequently, PROPACK is an efficient — both in terms of number of flops and storage 
requirement — and stable package for computing the dominant singular values and singular vectors 
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of a large sparse matrix. For information, the available documentation [36] reports a speedup 
factor of about ten over Matlab's 'svds'. Furthermore, the Fortran version of PRO PACK is about 
3-4 times faster than the Matlab version. Despite this significant speedup, we have only used the 
Matlab version but since the singular value shrinkage operator is by-and-large the dominant cost in 
the SVT algorithm, we expect that a Fortran implementation would run about 3 to 4 times faster. 

As for most SVD packages, though one can specify the number of singular values to compute, 
PRO PACK can not automatically compute only those singular values exceeding the threshold r. 
One must instead specify the number s of singular values ahead of time, and the software will 
compute the s largest singular values and corresponding singular vectors. To use this package, we 
must then determine the number Sk of singular values of Y^~^ to be computed at the /cth iteration. 
We use the following simple method. Let r^-i = rank{X^~^) be the number of nonzero singular 
values of X^~^ at the previous iteration. Set = r^-i + 1 and compute the first singular values 
of Y^~^. If some of the computed singular values are already smaller than r, then Sk is a right 
choice. Otherwise, increment by a predefined integer i repeatedly until some of the singular 
values fall below r. In the experiments, we choose i = 5. Another rule might be to repeatedly 
multiply Sk by a positive number — e.g. 2 — until our criterion is met. Incrementing by a fixed 
integer works very well in practice; in our experiments, we very rarely need more than one update. 

We note that it is not necessary to rerun the Lanczos iterations for the first vectors since they 
have been already computed; only a few new singular values {i of them) need to be numerically 
evaluated. This can be done by modifying the PROPACK routines. We have not yet modified 
PROPACK, however. Had we done so, our run times would be decreased. 

5.1.2 Step sizes 

There is a large literature on ways of selecting a step size but for simplicity, we shall use step sizes 
that are independent of the iteration count; that is Sk = ^ for k = 1,2,.... From Theorem 4.2, 
convergence for the completion problem is guaranteed (2.7) provided that < 6 < 2. This choice 
is, however, too conservative and the convergence is typically slow. In our experiments, we use 
instead 



i.e. 1.2 times the undersampling ratio. We give a heuristic justification below. 

Consider a fixed matrix A G R"i^"2^ Under the assumption that the column and row spaces of 
A are not well aligned with the vectors taken from the canonical basis of M"^ and M"^ respectively — 
the incoherence assumption in [13] — then with very large probability over the choices of 0, we have 



provided that the rank of A is not too large. The probability model is that is a set of sampled 
entries of cardinality m sampled uniformly at random so that all the choices are equally likely. In 
(5.2), we want to think of e as a small constant, e.g. smaller than 1/2. In other words, the 'energy' 
of A on 17 (the set of sampled entries) is just about proportional to the size of The near isometry 
(5.2) is a consequence of Theorem 4.1 in [13], and we omit the details. 



Now returning to the proof of Theorem 4.2, we see that a sufficient condition for the convergence 
of (2.7) is 



(1 - e)p WAfp < \\rn{A)\\l < (1 + e)p \\A\\l, p := ^/(nina), 



(5.2) 



3;^ > 



26\\X* - X'^Wl + d'^WVniX" - X'')\\l < -pWX" - X 
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compare (4.4), which is equivalent to 

Since ||'Pn(X)||ir < for any matrix X G M"i^"2^ jg gg^fg gelect 5 <2. But suppose that 

we could apply (5.2) to the matrix A = X* — X^. Then we could take 6 inversely proportional 
to p; e.g. with e = 1/4, we could take 6 < 1.6p~^. Below, we shall use the value 6 = 1.2p~^ which 
allows us to take large steps and still provides convergence, at least empirically. 

The reason why this is not a rigorous argument is that (5.2) cannot be applied to A = X* — X^ 
even though this matrix difference may obey the incoherence assumption. The issue here is that 
X* — X^ is not a fixed matrix, but rather depends on Q since the iterates {X''} are computed 
with the knowledge of the sampled set. 

5.1.3 Initial steps 

The SVT algorithm starts with = 0, and we want to choose a large r to make sure that the 
solution of (2.8) is close enough to a solution of (1.1). Define kQ as that integer obeying 

r 



G{ko-l,ko]. (5.3) 



S\\rn{M)\\2 

Since = 0, it is not difficult to see that 

X^ = 0, Y^ = k5Vn{M), k=l,...,ko. 

To save work, we may simply skip the computations of X^, . . . ,X^°, and start the iteration by 
computing X'^o+i from Y''° . 

This strategy is a special case of a kicking device introduced in [44]; the main idea of such 
a kicking scheme is that one can 'jump over' a few steps whenever possible. Just like in the 
aforementioned reference, we can develop similar kicking strategies here as well. Because in our 
numerical experiments the kicking is rarely triggered, we forgo the description of such strategies. 

5.1.4 Stopping criteria 

Here, we discuss stopping criteria for the sequence of SVT iterations (2.7), and present two possi- 
bilities. 

The first is motivated by the first-order optimality conditions or KKT conditions tailored to the 
minimization problem (2.8). By (2.14) and letting dYSoiY) = in (2.13), we see that the solution 
X* to (2.8) must also verify 

|^=^'<^'- (5.4) 
\P„{X - M) = 0, 

where 1^ is a matrix vanishing outside of Q'^. Therefore, to make sure that X^ is close to X*, it 
is sufficient to check how close (X'^, "K'^"^) is to obeying (5.4). By definition, the first equation in 
(5.4) is always true. Therefore, it is natural to stop (2.7) when the error in the second equation is 
below a specified tolerance. We suggest stopping the algorithm when 

\\Vn{X'^-M)y 
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where e is a fixed tolerance, e.g. 10~^. We provide a short heuristic argument justifying this choice 
below. 

In the matrix completion problem, we know that under suitable assumptions 

\\Vn{M)\\l>ip\\M\\l, 

which is just (5.2) applied to the fixed matrix M (the symbol x here means that there is a constant 
e as in (5.2)). Suppose we could also apply (5.2) to the matrix — M (which we rigorously cannot 
since X^ depends on 0), then we would have 

WTniX"" -M)\\l-:<p\\X^ -M\\l, (5.6) 

and thus 

\\Vn{X^-M)\\F ^ \\X^-M\\f 
\\Vn{M)\\F \\M\\f ■ 

In words, one would control the relative reconstruction error by controlling the relative error on 
the set of sampled locations. 

A second stopping criterion comes from duality theory. Firstly, the iterates X^ are generally 
not feasible for (2.8) although they become asymptotically feasible. One can construct a feasible 
point from X^ by projecting it onto the affine space {X : Vn{X) = Vn{M)} as follows: 

X'' = X'' + Vn{M - X''). 

As usual let fri^) = ''"ll^ll* + ^II-^IIf and denote by p* the optimal value of (2.8). Since X'^ is 
feasible, we have 

Secondly, using the notations of Section 2.4, duality theory gives that 

ak := <7o(l^'~') = /:(X^l^'=-l) < p\ 

Therefore, bk — ak is an upper bound on the duality gap and one can stop the algorithm when this 
quantity falls below a given tolerance. 

For very large problems in which one holds X^ in reduced SVD form, one may not want to 
compute the projection X'^ since this matrix would not have low rank and would require signifi- 
cant storage space (presumably, one would not want to spend much time computing this projection 
either). Hence, the second method only makes practical sense when the dimensions are not pro- 
hibitively large, or when the iterates do not have low rank. 

5.1.5 Algorithm 

We conclude this section by summarizing the implementation details and give the SVT algorithm 
for matrix completion below (Algorithm 1). Of course, one would obtain a very similar structure 
for the more general problems of the form (3.1) and (3.4) with linear inequality constraints. For 
convenience, define for each nonnegative integer s < min{?ii, 722}, 

[^7^5]^y^^, k = i,2,..., 

where = [ui, . . . , Ug] and = [vf, . . . , i?^] are the first s singular vectors of the matrix 1^^, 
and 5]*^ is a diagonal matrix with the first s singular values , . . . , cr^ on the diagonal. 
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Algorithm 1: Singular Value Thresholding (SVT) Algorithm 



Input: sampled set and sampled entries Vq{M), step size 5, tolerance e, parameter 
r, increment i, and maximum iteration count k^g^x 
Output: X°P* 

Description: Recover a low-rank matrix M from a subset of sampled entries 

1 Set = ko6Vn{M) {Uq is defined in (5.3)) 

2 Set ro = 

3 for A; = 1 to A;jnax 

4 Set Sk = Tfc-i + 1 

5 repeat 

6 Compute [C/^-\S^-i,F'=^i],, 

7 Set Sfc = Sfc + £ 

8 until (T^s~\ < T 

9 Set Tk = max{j : > t} 

10 Set X'' = y." , (cj^-i - r)u^^ii7'=-i 

11 



12 



if llPnCX'^ - M)||F/||PnM||F < e then break 



Set 

13 end /or /c 

14 Set X°P* = X 



_ )0 if(i,j)0S^, 
\y^^-i + 5(M,,-X,^.) if(z,j)Gfi 



k 



5.2 Numerical results 

5.2.1 Linear equality constraints 

Our implementation is in Matlab and all the computational results we are about to report were 
obtained on a desktop computer with a 1.86 GHz CPU (dual core with Matlab's multithreading 
option enabled) and 3 GB of memory. In our simulations, we generate n x n matrices of rank r 
by sampling two n x r factors Ml and Mpi independently, each having i.i.d. Gaussian entries, and 
setting M = M^M^I^ as it is suggested in [13]. The set of observed entries Q is sampled uniformly 
at random among all sets of cardinality m. 

The recovery is performed via the SVT algorithm (Algorithm 1), and we use 

WVniX'' - M)\\F/\\VnM\\F < lO'^ (5.7) 

as a stopping criterion. As discussed earlier, the step sizes are constant and we set 6 = 1.2p~^. 
Throughout this section, we denote the output of the SVT algorithm by X°p*. The parameter r 
is chosen empirically and set to r = 5n. A heuristic argument is as follows. Clearly, we would like 
the term T||iW"||^, to dominate the other, namely, ^||A^"|||,. For products of Gaussian matrices as 
above, standard random matrix theory asserts that the Frobenius norm of M concentrates around 
n^/r, and that the nuclear norm concentrates around about nr (this should be clear in the simple 
case where r = 1 and is generally valid). The value t = 5n makes sure that on the average, the 
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Unknown M 


Computational results 


size (n x n) 


rank (r) 


m /dr 




time(s) 


# iters 


relative error 




10 


6 


0.12 


23 


117 


1.64 X 10~^ 


1,000 X 1,000 


50 


4 


0.39 


196 


114 


1.59 X 10"^ 




100 


3 


0.57 


501 


129 


1.68 X 10"^ 




10 


6 


0.024 


147 




1.16 X iU 


5, 000 X 5, 000 


50 


5 


0.10 


950 


108 


1.61 X 10"^ 




100 


4 


0.158 


3,339 


123 


1.72 X 10-^ 




10 


6 


0.012 


281 


123 


1.73 X 10""^ 


10,000 X 10,000 


50 


5 


0.050 


2,096 


110 


1.65 X 10-^ 




100 


4 


0.080 


7,059 


127 


1.79 X 10"^ 


20, 000 X 20, 000 


10 


6 


0.006 


588 


124 


1.73 X 10-4 






50 


5 


0.025 


4,581 


111 


1.66 X 10-4 


30, 000 X 30, 000 


10 


6 


0.004 


1,030 


125 


1.73 X 10-4 



Table 1: Experimental results for matrix completion. The rank r is the rank of the unknown 
matrix M, m/dr is the ratio between the number of sampled entries and the number of 
degrees of freedom in an n x n matrix of rank r (oversampling ratio), and mjv? is the fraction 
of observed entries. All the computational results on the right are averaged over five runs. 

value of r||iVi"||* is about 10 times that of as long as the rank is bounded away from the 

dimension n. 

Our computational results are displayed in Table 1. There, we report the run time in seconds, the 
number of iterations it takes to reach convergence (5.7), and the relative error of the reconstruction 

relative error = ||X°p* - M||F/||M||i;', (5.8) 

where iW is the real unknown matrix. All of these quantities are averaged over five runs. The table 
also gives the percentage of entries that are observed, namely, mjv? together with a quantity that 
we may want to think as the information oversampling ratio. Recall that an n x n matrix of rank 
r depends upon ■= r[2n — r) degrees of freedom. Then m/dr is the ratio between the number of 
sampled entries and the 'true dimensionality' of an n x n matrix of rank r. 

The first observation is that the SVT algorithm performs extremely well in these experiments. 
In all of our experiments, it takes fewer than 200 SVT iterations to reach convergence. As a 
consequence, the run times are short. As indicated in the table, we note that one recovers a 
1, 000 X 1, 000 matrix of rank 10 in less than a minute. The algorithm also recovers 30, 000 x 30, 000 
matrices of rank 10 from about 0.4% of their sampled entries in just about 17 minutes. In addition, 
higher-rank matrices are also efficiently completed: for example, it takes between one and two 
hours to recover 10, 000 x 10, 000 matrices of rank 100 and 20, 000 x 20, 000 matrices of rank 50. We 
would like to stress that these numbers were obtained on a modest CPU (1.86GHz). Furthermore, 
a Fortran implementation is likely to cut down on these numbers by a multiplicative factor typically 
between three and four. 

We also check the validity of the stopping criterion (5.7) by inspecting the relative error defined 
in (5.8). The table shows that the heuristic and nonrigorous analysis of Section 5.1 holds in practice 
since the relative reconstruction error is of the same order as ||'Pn(X°P* — iVf)||i?/||'Pt7iW"||i? ~ 10-^. 
Indeed, the overall relative errors reported in Table 1 are all less than 2 x 10"^. 
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We emphasized all along an important feature of the SVT algorithm, which is that the matrices 
have low rank. We demonstrate this fact empirically in Figure 1, which plots the rank of 
versus the iteration count k, and does this for unknown matrices of size 5, 000 x 5, 000 with 
different ranks. The plots reveal an interesting phenomenon: in our experiments, the rank of X^ 
is nondecreasing so that the maximum rank is reached in the final steps of the algorithm. In fact, 
the rank of the iterates quickly reaches the value r of the true rank. After these few initial steps, 
the SVT iterations search for that matrix with rank r minimizing the objective functional. As 
mentioned earlier, the low-rank property is crucial for making the algorithm run fast. 



10 



50 



100 



Figure 1: Rank of X as a function k when the unknown matrix M is of size 5, 000 x 5, 000 
and of rank r. 

Finally, we demonstrate the results of the SVT algorithm for matrix completion from noisy 
sampled entries. Suppose we observe data from the model 

Bij=Mij + Zij, (5.9) 

where Z is a zero-mean Gaussian white noise with standard deviation a. We run the SVT algorithm 
but stop early, as soon as X^ is consistent with the data and obeys 

\\Vi,{X''-B)\\l<{l + e)ma\ (5.10) 

where e is a small parameter. Our reconstruction M is the first X^ obeying (5.10). The results 
are shown in Table 2 (the quantities are averages of 5 runs). Define the noise ratio as 

\\VniZ)\\F/\\Vnm\\F, 

and the relative error by (5.8). From Table 2, we see that the SVT algorithm works well as the 
relative error between the recovered and the true data matrix is just about equal to the noise ratio. 

The theory of low-rank matrix recovery from noisy data is nonexistent at the moment, and is 
obviously beyond the scope of this paper. Having said this, we would like to conclude this section 
with an intuitive and nonrigorous discussion, which may explain why the observed recovery error 
is within the noise level. Suppose again that M obeys (5.6), namely, 

\\Vn{M - M)\\l^p\\M - M\\l. (5.11) 
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noise ratio 


Unknown matrix M 


Computational results 




size (n x n) 


rank (r) 


m/dr 


m In^ 

1 


time(s) 


# iters 


relative error 






10 


6 


0.12 


10.8 


51 


0.78 X 10"^ 




1,000 X 1,000 


50 


4 


0.39 


87.7 


48 


0.95 X 10-2 






100 


3 


0.57 


216 


50 


1.13 X 10-2 






10 


6 


0.12 


4.0 


19 


0.72 X 10-1 


10-1 


1,000 X 1,000 


50 


4 


0.39 


33.2 


17 


0.89 X 10-1 






100 


3 


0.57 


85.2 


17 


1.01 X 10-1 






10 


6 


0.12 


0.9 


3 


0.52 


1 


1,000 X 1,000 


50 


4 


0.39 


7.8 


3 


0.63 






100 


3 


0.57 


34.8 


3 


0.69 



Table 2: Simulation results for noisy data. The computational results are averaged over five 
runs. 

As mentioned earlier, one condition for this to happen is that M and NL have low rank. This is 
the reason why it is important to stop the algorithm early as we hope to obtain a solution which 
is both consistent with the data and has low rank (the limit of the SVT iterations, lim^^oo X'^, 
will not generally have low rank since there may be no low-rank matrix matching the noisy data). 
From 

||Pq(M - M)\\f < \\Vn{M - B)\\f + \\Vn{B - M)\\f, 

and the fact that both terms on the right-hand side are on the order of V ma"^, we would have 
p||Ai" — -/Vf|||n = 0{ma'^) by (5.11). In particular, this would give that the relative reconstruction 
error is on the order of the noise ratio since ||'Pn(-M")|||;, x p||iW"||^ — as observed experimentally. 

5.2.2 Linear inequality constraints 

We now examine the speed at which one can solve similar problems with linear inequality constraints 
instead of linear equality constraints. We assume the model (5.9), where the matrix M of rank 
r is sampled as before, and solve the problem (3.8) by using (3.10). We formulate the inequality 
constraints in (3.8) with Eij = cj so that one searches for a solution M with minimum nuclear 
norm among all those matrices whose sampled entries deviate from the observed ones by at most 
the noise level cr.^ In this experiment, we adjust a to be one tenth of a typical absolute entry of 
M, i.e. a = 0.1 X^j^gfj \Mij\/m, and the noise ratio as defined earlier is 0.780. We set n = 1,000, 
r = 10, and the number m of sampled entries is five times the number of degrees of freedom, 
i.e. m = hdr- Just as before, we set r = 5n, and choose a constant step size 6 = 1.2p~^ . 

The results, reported in Figure 2, show that the algorithm behaves just as well with linear 
inequality constraints. To make this point, we compare our results with those obtained from 
noiseless data (same unknown matrix and sampled locations). In the noiseless case, it takes about 
150 iterations to reach the tolerance e = 10"^ whereas in the noisy case, convergence occurs in 
about 200 iterations (Figure 2(a)). In addition, just as in the noiseless problem, the rank of the 
iterates is nondecreasing and quickly reaches the true value r of the rank of the unknown matrix 

^This may not be conservative enough from a statistical viewpoint but this works well in this case, and our 
emphasis here is on computational rather than statistical issues. 
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(a) 



(b) 



(c) 



Figure 2: Computational results of the algorithm applied to noisy (linear inequality con- 
straints as in (3.8)) and noiseless data (equality constraints). The blue (resp. red) color is 
used for the noisy (resp. noiseless) experiment, (a) Plot of the reconstruction errors from 
noisy and noiseless data as a function of the iteration count. The thin line is the residual 
relative error \\Vn{X^ — M)\\p /\\Vq,{M)\\f and the thick line is the overall relative error 
\\X^ — M\\f /\\M\\p. (b) Rank of the iterates as a function of the iteration count, (c) Time it 
takes to compute the singular value thresholding operation as a function of the iteration count. 
The computer here is a single-core 3.00GHz Pentium 4 running Matlab 7.2.0. 

M we wish to recover (Figure 2(b)). As a consequence the SVT iterations take about the same 
amount of time as in the noiseless case (Figure 2(c)) so that the total running time of the algorithm 
does not appear to be substantially different from that in the noiseless case. 

We close by pointing out that from a statistical point of view, the recovery of the matrix M 
from undersampled and noisy entries by the matrix equivalent of the Dantzig selector appears to 
be accurate since the relative error obeys \\M — M\\f /\\M\\f = 0.0769 (recall that the noise ratio 
is about 0.08). 

6 Discussion 



This paper introduced a novel algorithm, namely, the singular value thresholding algorithm for 
matrix completion and related nuclear norm minimization problems. This algorithm is easy to 
implement and surprisingly effective both in terms of computational cost and storage requirement 
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when the minimum nuclear-norm solution is also the lowest-rank solution. We would like to close 
this paper by discussing a few open problems and research directions related to this work. 

Our algorithm exploits the fact that the sequence of iterates {-X^'^} have low rank when the 
minimum nuclear solution has low rank. An interesting question is whether one can prove (or 
disprove) that in a majority of the cases, this is indeed the case. 

It would be interesting to explore other ways of computing T>t-(Y) — in words, the action of 
the singular value shrinkage operator. Our approach uses the Lanczos bidiagonalization algorithm 
with partial reorthogonalization which takes advantages of sparse inputs but other approaches are 
possible. We mention two of them. 

1. A series of papers have proposed the use of randomized procedures for the approximation 
of a matrix Y with a matrix Z of rank r [38,41]. When this approximation consists of the 
truncated SVD retaining the part of the expansion corresponding to singular values greater 
than r, this can be used to evaluate T>r(Y). Some of these algorithms are efficient when the 
input Y is sparse [41], and it would be interesting to know whether these methods are fast 
and accurate enough to be used in the SVT iteration (2.7). 

2. A wide range of iterative methods for computing matrix functions of the general form f{Y) 
are available today, see [34] for a survey. A valuable research direction is to investigate 
whether some of these iterative methods, or other to be developed, would provide powerful 
ways for computing DriY). 

In practice, one would like to solve (2.8) for large values of r. However, a larger value of r 
generally means a slower rate of convergence. A good strategy might be to start with a value of 
r, which is large enough so that (2.8) admits a low-rank solution, and at the same time for which 
the algorithm converges rapidly. One could then use a continuation method as in [50] to increase 
the value of r sequentially according to a schedule ro,ri, . . ., and use the solution to the previous 
problem with r = rj_i as an initial guess for the solution to the current problem with r = Tj (warm 
starting). We hope to report on this in a separate paper. 
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